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We present experimental results on the velocity statistics of a uniformly heated granular fluid, 
in a quasi-2D configuration. We find the base state, as measured by the single particle velocity 
distribution P(c), to be universal over a wide range of filling fractions and only weakly dependent 
on all other system parameters. There is a consistent overpopulation in the distribution's tails, 
which scale as P oc exp(— A x c~^^'^). More generally, P{c) deviates from a Maxwell-Boltzmann 
by a second order Sonine polynomial with a single adjustable parameter, in agreement with recent 
theoretical analysis of inelastic hard spheres driven by a stochastic thermostat. To our knowledge, 
this is the first time that Sonine deviations have been measured in an experimental system. 



I. INTRODUCTION 

The study of granular flows has received a recent up- 
surge of interest in physics Jj. This has been motivated 
by both the relevance of such flows to a wide range of in- 
dustrial and geological processes, and by the realization 
that granular materials provide an excellent test bed for a 
number of fundamental question in the context of modern 
fluid dynamics and non-equilibrium statistical mechanics 
Q . Granular media are ensembles of macroscopic parti- 
cles, in which kinetic energy {Ke) can be preferentially 
injected in the N ~ 10 — 10^" translational and rotational 
degrees of freedom of the center of mass of the particles 
without effecting the Na — 10^^ internal degrees of free- 
dom of each particle. Therefore the granular temperature 
T = Ke/N is 10 — 20 orders of magnitude larger than 
the molecular temperature r = Ke/N a- In this sense the 
molecular temperature is irrelevent. Moreover, granular 
materials are intrinsically dissipative since energy from 
the Ke of the particles is lost (converted) through in- 
elastic collisions and frictional contacts to the molecular 
temperature. Hence, any dynamical study of a granular 
media requires an energy input, which typically takes the 
form of vibration or shear [3|, This interplay between 
energy injection and dissipation forces the system to be 
far from equilibrium and granular fluids are therefore not 
expected to behave similarly to the equilibrium counter- 
parts of classical fluids 0]. 

However, we have recently shown that the structure 
of a uniformly heated quasi-2D granular fluid is iden- 
tical to that found in simulations of equilibrium hard 
disksf^. Further, we find that the mean-squared dis- 
placement of grains shows diffusive behavior at low den- 
sity and caging behavior as density is increased until 
particles are completely confined above freezing point 
just as in ordinary molecular fluids^. These surpris- 
ing results provide a partial experimental justification 
for the theoretical treatment of granular flow by anal- 



ogy with molecular fluids using standard statistical me- 
chanics and kinetic theory_^ A large number of kinetic 
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221 . I23l |2j| have been pro- 
duced starting in 1980's. The results of these theories are 
balance equation analogous to the Navier-Stokes equa- 
tions for the density, n, the average velocity, Vq, and the 
average kinetic energy per particle or temperature, T, 
as well as, equations of state and transport coefhcients, 
which are functions of n,vo, and T. The heart of these 
formulations is the single particle distribution function 
/(I) (r,v,<), which measures the probability of a particle 
a position r having velocity v. The field equations are 
derived from moments of /^^^ (r,v, i). For example, the 
density, n, is the 0th moment, the average velocity, Vq, 
is the 1st moment, and the temperature, T, is the 2nd 
moment. The transport coefficients are deterimincd from 
non-equilibrium perturbations of f^^^ (r,v,t). Generally 
the equilibrium or steady-state distribution is assumed 
to be indcpent of space and to have Maxwell-Boltzmann 

form /(^^ = = Aexp 2t"^ ) ' a-ltbough, recent 
work has explored the use of other base states or steady- 
state distributions [111,!!!!!!]. 

Two important questions one may ask are: 1) What 
is the base state /^^^ for a driven granular fluid? and 2) 
Does it have a universal form?. If it is characterized by a 
small number of parameters such as the moments of the 
f^^^ distribution as in the case for ideal fluids, one would 
be able to advance a great deal in developing predictive 
continuum models for granular flows since the theoretical 
machinery from kinetic theory could be readily applied. 
The careful experimental study and characterization of 
this base state is therefore important in order to develop 
a theory for granular media similar to that of regular 
molecular liquids. 

There has been a substantial amount of experimen- 
tal, numerical, and theoretical work to address these is- 
sues in configurations where the energy input perfectly 
balances the dissipation such that the system reaches a 
steady state, albeit far from equilibrium. These non- 
equilibrium steady states are simplified realizations of 
granular flows more amenable to analysis, but the in- 
sight gained from their study should be helpful in the 
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tackling of other more complex scenarios. One fea- 
ture that has been consistantly established in experi- 
ments is that single particle velocity distribution func- 
tions deviate from the Maxwell-Boltzmann distribution 
dS m, [13, m, S H iii- In particular, the tails (i.e. 
large velocities) of the experimental distributions exhibit 
a considerable overpopulation and have been shown to 
scale as f'-^^^Vi) ~ exp[{vi/y/Ti)^^'^] where Vi is a velocity 
component and Tj = (vf) is the granular temperature. 
This behavior is in good agreement with both numerical 
[3^ and theoretical [s^ predictions. Even though the 
tails correspond to events with extremely low probabili- 
ties, they increase the variance of the distribution. The 
variance of the distribution is T, the average kinetic en- 
ergy of the particles. Using a Gaussian to represent this 
type of distribution leads to major discrepancy in the 
region of high probability at the central part of the dis- 
tribution (low velocities) [3l| which have, thus far, been 
greatly overlooked in experimental work. Note that re- 
quiring the variance of f^^\vi) to be the granular tem- 
perature greatly constrains the functional form that the 
distribution can take. 

A model system that has been introduced to study 
this question is the case of a homogeneous granular gas 
heated by a stochastic thermostat, i.e. an ensemble 
of inelastic particles randomly driven by a white noise 
energy source [sj]. Recently, there have been many 
theoretical and numerical studies on this model system 
[H,!!!,!!!,!!!, [13, nil where the steady state velocity dis- 
tribution have been found to deviate from the Maxwell- 
Boltzmann distribution, van Noije and Ernst [s^ stud- 
ied these velocity distributions based on approximate so- 
lutions to the inelastic hard sphere Enskog-Boltzmann 
equation by an expansion in Sonine polynomials. The 
results of their theoretical analysis has been validated 
by numerical studies using both molecular dynamics 133 
simulations and direct simulation Monte Carlo [H, |38| |. 
The use Sonine corrections is particular attractive since 
it leaves the variance of the resulting velocity distribution 
unchanged but leads to a non-Gaussian fourth moment 
or kurtosis K 3. 

We have addressed the above issues in a well con- 
trolled experiment in which we are able to perform pre- 
cision measurements of the velocity distributions of a 
uniformly heated granular fluid. A novel feature of our 
experimental technique is that we are able to generate 
macroscopic random walkers over a wide range of filling 
fractions and thermalize the granular particles in a way 
analogous to a stochastic thermostat. This is an ideal 
system to test the applicability of some of the kinetic 
theory ideas mentioned above. In our earlier study Q, 
we focused on the structural configurations of this granu- 
lar fluid and revealed striking similarity to those adopted 
by a fluid of equilibrium hard disks. Here we concentrate 
on the dynamical aspects of this same experimental sys- 
tem, as measured by the single particle velocity distri- 
bution, P{c), and observe a marked departure from the 
equilibrium behavior, i.e. from the Maxwell-Botlzmann 



distribution. We quantify these deviations from equilib- 
rium and show that they closely follow the predictions, 
in the kinetic theory framework, of the theoretical and 
numerical work mentioned above. In particular, we find 
a consistent overpopulation in the distribution's high en- 
ergy tails, which scale as stretched exponentials with ex- 
ponent —3/2. In addition, we experimentally determine 
the deviations from a Gaussian of the full distributions 
using a Sonine expansion method and find them to be 
well described by a second order Sonine polynomial cor- 
rection. We establish that the entire experimental veloc- 
ity distribution functions are universal and well described 
by 

Pic) = fMB (1 + a2(l/2c4 - 3/22 ^ 3/8)^ 

where 02 is the non-gaussianity parameter that can be 
related to the Kurtosis of the distribution and the forth- 
order polynomial is the Sonine polynomial of order 2. 
This functional form of P(c) was found to be valid over a 
wide range of the system parameters. To our knowledge, 
this is the first time that the Sonine corrections of the 
velocity distributions are measured in an experimental 
system, in agreement with analytical predictions. 

This paper is organized as follows. In Sec. [TT]we briefly 
review the theoretical framework of kinetic theory along 
with the steady state solution to the Enskog-Boltzmann 
equation using the Sonine polynomial expansion method. 
In Sec. mil we present our experimental apparatus and 
describe how the homogeneous granular gas with random 
heating is generated. In Sec. llVl we study the thermaliza- 
tion mechanism by analyzing the trajectories of a single 
particle in the granular cell. In Sec.,|Vl using the concept 
of granular temperature, we then quantify the dynamics 
of the experimentally obtained non-equilibrium steady 
states, for a number of parameters, namely: the filling 
fraction, the driving parameters (frequency and acceler- 
ation) and the vertical gap of the cell. In Sec. I VII we turn 
to a detailed investigation of the Probability Distribution 
Functions of velocities as a function of the system pa- 
rameters. In particular, we quantify the deviations from 
Maxwell-Boltzmann behavior at large velocities (tails of 
the distributions) - Sec. NTM In Sec. \VTB\ we extend 
the deviation analysis to the full range of the distribution 
using an expansion method that highlights deviations at 
low velocities and allows us to make a connection with 
Sonine polynomials. We finish in Sec. IVIII with some 
concluding remarks. 

II. BRIEF REVIEW OF THEORY 

We briefiy review the key features of the kinetic theory 
for granular media pertinent to our study. The number 
of particles in a volume element, dr, and velocity ele- 
ment, dv, centered at position r and velocity v is given 
by /^(r, V, i)drdv, where f^{r,v,t) is the single particle 
distribution function. Continuum fields are given as av- 
erages over /^(r,v,t). For instance, the number density. 
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n, average velocity, Wo, and granular temperature, T, are 
given respectively by. 



n(r, t) 



,t)dv. 



Vo(r,i) = /i(r,v,i)vdv, 

r(r,i) = ^ / /i(r,v,t)(v-v„)2dv, 
nd J 



(2) 



(3) 



(4) 



where d is the number of dimensions. It is important 
to note that the granular temperature T is not the ther- 
modynamic temperature but, by analogy, the kinetic en- 
ergy per macroscopic particles (explained in detail in Sec. 
W} . For the spatially homogeneous and isotropic case, 
we refer to the single particle distribution function as 
f^{v,t) and consider only a single component of the ve- 
locity, i.e. f^^\vi,t). From now on we drop the index 
i. It is also convenient to introduce a scaled distribution 
function f^{c,t) where the velocity is scaled by a char- 
acteristic velocity such that, c = v/ yj2T{t), where T{t) 
is the granular temperature and equal to the variance of 

The stocastically heated single particle distribution 
function f{c,t) satisfies the Enskog-Boltzmann equation. 



9/ 
dt 



QifJ)+^Fpf, 



(5) 



where, Q, is the collision operator, which accounts for the 
inelastic particle interactions and the Fokker-Plank op- 
erator, J^FP, accounts for the stochastic forcing [ssHsoj. 
We are interested in a stationary solution of Eqn. ^ 
where the heating exactly balances the loss of energy due 
to collisions and the temperature becomes time indepen- 
dent, van Noije and Ernst [s^ obtained steady state 
solutions to Eqn. ([5|) by taking the series expansion of 
/(c) away from a Maxwell-Boltzmann, /mb, i-e- 



/(c) = /ms(c) 




.5p(c2) 



(6) 



p=i 



in terms of the Soninc polynomials Sp{c^) where Op is 
a numerical coefficient for the p'^'-order. The deviations 
from the MB distribution are thus expressed in terms of 
an expansion on Sonine polynomials. The Sonine poly- 
nomials, which are knows as associated Laguerre polyno- 
mials in different contexts, are defined as, 

_ ' ^l)"(p+l/2)! 

These are orthogonal accordingly to the relationship, 

r c^e-^'s,{^)Sm{^)dc^U,J^^. (8) 
.In ^ n- 



To further simply the calculation, terms higher than 0(2) 
are typically neglected, such that. 



where, 



/(c)-/A/B(c){l + a252(c2)} 



/a/b(c) = — =ea;p(-c 
a/tt 



is the Maxwell-Boltzmann distribution and 
a I i\ 1 4 3 2 3 

^2(C) = -C --C 



(9) 



(10) 



(11) 



is the 0(2) Sonine polynomial for 2D. The first Sonine 
coefficient, a\ , is zero according to the definition of tem- 
perature [40|. The second Sonine coefficient, is the 
first non-trivial Sonine coefficient and hence the first non- 
vanishing correction to the MB distribution and can be 
related to the Kurtosis K = 3 (02 + 1) of the distribution. 



III. OUR EXPERIMENTS 

In Fig. [1] we present a schematic diagram of our ex- 
perimental apparatus, which is adapted from a design 
introduced by Olafson and Urbach [2^, [2^ . An ensem- 
ble of stainless steel spheres (diameter D = 1.19mm) 
were vertically vibrated in a shallow cylindrical cell, at 
constant volume conditions. The cell consisted of two 
parallel glass plates set horizontally and separated by a 
stainless steel annulus. This annulus had an inner di- 
ameter of 101.6mm and, unless otherwise stated, was 
1.905mm thick, which set the height of the cell to be 1.6 
particle diameters. This constrained the system to be 
quasi-two dimensional such that each sphere could not 
fully overlap. However the 2D projection of the particles 
have a minimum separation of r„iin — O-SD (rather than 
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FIG. 1: Schematic diagram of the experimental apparatus. 
The ensemble of stainless steel spheres was sandwiched in be- 
tween two glass plates, separated by a 1.905mm thick annulus 
(dark grey). Top glass plate was always optically flat. Bot- 
tom glass plate could either be (a) optically flat or (b) rough 
by sand-blasting. 
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Tmin = D for a strictly two-dimensfonal system). Hence, 
the extreme situation of a maximum overlap between two 
particles {S — D — Tmin — 0.2D) occurred when one of 
the particles of the adjacent pair was in contact with the 
bottom plate and the other in contact with the top plate, 
as shown in Fig. [Ija). We have also performed experi- 
ments to explore the effect of changing the height of the 
cell. 

The top disk of the cell was made out of an optically 
flat borofloat glass treated with a conducting coating of 
ITO to eliminate electrostatic effects. For the bottom 
disk of the cell, two variants were considered: a flat glass 
plate identical to the top disk - Fig. [Ija) - and a rough 
borosilicate glass plate - Fig. [TJb). In the second case, 
the plate was roughened by sand-blasting, generating 
random structures with lengthscales in the range of 50- 
500/^m. The use of a rough bottom glass plate improves 
on the setup of Olafsen and Urbach 26] , who used a flat 
plate. As it will be shown in Sec. IIVI and Appendix \^ 
it had the considerable advantage of effectively random- 
izing the trajectories of single particles under vibration, 
allowing a considerably wider range of filling fractions to 
be explored. 

The horizontal experimental cell was vertically vi- 
brated, sinusoidally, via an electromagnetic shaker 
(VGlOO-6 Vibration Test System). The connection of the 
shaker to the cell was done via a robust rectangular linear 
air-bearing which constrained the motion to be unidirec- 
tional. The air-bearing ensured filtering of undesirable 
harmonics and non-vertical vibrations due to its high 
stiffness provided by high-pressure air flow in between 
the bearing surfaces. Moreover, the coupling between 
the air-bearing and the shaker consisted of a thin brass 
rod (25.4mm long and 1.6mm diameter). This rod could 
slightly flex to correct any misalignment present in the 
shaker/bearing system, while being sufficiently rigid in 
the vertical direction to fully transmit the motion. This 
ensemble - shaker/brass rod/air-bearing - ensured a high 
precision vertical oscillatory driving. 

The forcing parameters of the system were the fre- 
quency, /, and amplitude. A, of the sinusoidal oscilla- 
tions. From these, it is common practice to construct a 
non-dimensional acceleration parameter. 



ter was the height, h, of the experimental cell which we 
varied from 1.3D <h < 2.3D. 

The granular cell was set horizontal in order to min- 
imize compaction effects, inhomogeneities and density 
gradients which otherwise would be induced by gravity. 
This way, a wide range of filling fractions, < 4> < 0.8, 
could be accurately explored by varying the number of 
spheres in the cell, N, down to a resolution of single 
particle increments. Moreover, as we shall show in Sec. 
lIVfTVl we were able to attain spatially uniform driving of 
the spheres in the cell due to the use of the rough glass 
bottom plate. 

The dynamics of the system was imaged from above 
by digital photography using a grayscale DALSA CA-D6 
fast camera, at 840 frames per second. The granular layer 
was illuminated from below, in a transmission configura- 
tion, by a 16 X 16 array of high intensity LEDs. In this 
arrangement the particles obstruct the light source and 
appear as dark circles in a bright background. 

We have developed particle tracking software based on 
a two-dimensional least squares minimization of the indi- 
vidual particle positions, which is able to resolve position 
to sub-pixel accuracy. By focusing on a (15 x 15)mm^ 
imaging window located in the central region of the full 
cell we were typically able to achieve resolutions of 1 /20- 
1/10 of a pixel (which corresponds to 2.5~5/im). From 
the trajectories of the particles we could easily calculate 
an approximately instantaneous discretized velocity as, 

\v\t, +At)\^ ^{vAt, + At)r + vy{{t, + At))2 (14) 



with 



x'{tj + At) - x'it) 



At 



and v„ 



r 
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(12) 



j'{tj+At)-y'{t) 
At 
. 

for the i*'' particle in an experimental frame at time tj 
where At ~ 1.19/is is the time interval between two 
frames and (x, y) are the cartesian coordinates on the 
horizontal plane. For the remainder of this paper we 
shall be interested in the nature of the probability dis- 
tribution function of velocities, P{v), of the driven parti- 
cles. Statistics for calculating the P{v) distributions were 
constructed by analyzing 2048 frames at an acquisition 
rate of 840 frames per second which corresponded to an 
acquisition real time of 2.4381s. 



where g is the gravitational acceleration. We worked 
within the experimental ranges (10 < / < 100)Hz and 
1 < r < 6. The first control parameter was the filling 
fraction of the steel spheres in the cell defined as, 



NTT{D/2f 



(13) 



where N is the number of spheres (with diameter D ~ 
1.19mm) in the cell of radius R — 50.8mm. The filling 
fraction 4> is therefore defined in the projection of the cell 
onto the horizontal plane. The second control parame- 



IV. SINGLE PARTICLE DRIVING 

We first concentrate on the case of the driving of a 
single particle in the vibrating experimental cell. For 
each experimental run, a single sphere was positioned in 
the field of view of the camera with the aid of a strong 
magnet temporarily placed beneath the imaging window. 
The magnet was then removed and the trajectory of the 
sphere immediately recorded. This was repeated 100 
times to obtain statistics. Typical trajectories for two 
runs with a flat and rough bottom glass plates are shown 
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in Fig. [UJa) and (b), respectively. In both cases, the cell 
was vibrated at a frequency of / = 50Hz and dimen- 
sionless acceleration of F = 4.0. During each cycle of the 
driving, the motion of the single particle was thermalized 
due to collisions with both top and bottom glass plates. 

By comparing Fig. [2ta) with Fig. [2tb), it is clear 
that the particle excursions obtained by using a rough 
bottom plate are considerable larger than those for a flat 
plate. This can be readily quantified by calculating the 
Probability Density Functions of the velocities of single 
particles, P{v), which we plot in Fig. [U^c), for both cases. 
We now define the average kinetic energy, also known as 
granular temperature^ of a single particle as the variance 
of its velocity distribution. 



(16) 



where and Vy are the two orthogonal components of 
the velocity in the 2D horizontal plane and the brack- 
ets (.) denote time averages for the timeseries of the 
velocity components for 100 particle trajectories. Note 
that since we are dealing with monodisperse spheres, 
in this definition of kinetic temperature, the mass of 
the particle is taken to be unity. We have measured 
9.53 X 10~'*m^s~^ if a rough bottom glass plate 
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FIG. 2: Typical trajectories for single particle driven using a 
(a) rough or (b) flat bottom glass plate. Each of the trajec- 
tories is 2.438s long, (c) Probability distribution function of 
velocities for a single particles using rough (o) and flat (□) 
bottom glass plate. Driving parameters: / = hQHz, V = 4.0. 



is used and Tji^^ = 2.51 x 10~^m^s~^ for the case of the 
fiat bottom glass plate. This yields a temperature ratio 
between the two cases of T^oughf^fiat — 

Indeed, this is a signature that the structures in the 
rough plate were considerably more effective than the flat 
plate in transferring and randomizing the momentum of 
the steel spheres in the vertical direction (due to the sinu- 
soidal vertical oscillation of the cell) into the horizontal 
plane. The velocity of the steel sphere was randomized 
each time it collided with the peaks and valleys randomly 
distributed across the the sand-blasted rough glass sur- 
face. We can therefore regard this system as a spatially 
uniform heater of the granular particles, thereby gener- 
ating macroscopic random walkers. For this reason, and 
as it will become more clear in the next Sec, where we 
will further comment on the issue of isotropy, we shall 
focus our experimental study in using the rough bottom 
glass plate. 



V. DRIVEN MONOLAYERS: GRANULAR 
TEMPERATURE 

Having investigated the dynamics of single particles in 
the granular cell, we now turn to the study of granu- 
lar monolayers at higher filling fractions. In Fig. [3] we 
present typical configurations, at three values of 0, for 
a granular layer driven at / = 50Hz and F 4 with 
a rough bottom plate. The snap-shot in Fig. OJa), for 
(j) = 0.34, corresponds to a dilute state in which the par- 
ticles perform large excursions in between collisions as 
they randomly diffuse across the cell. If the filling frac- 
tion is increased, as in the example of the frame shown 
in Fig. mb) for cj) = 0.67, one observes a higher collision 
rate characteristic of a dense gaseous regime. For even 
higher values of filling fraction the spheres ordered into 
an hexagonally packed arrangement and became locked 
into the cage formed by its six neighbors. The system is 
then said to be crystalized as shown in the typical frame 
presented in Fig. [3]Jc) for cj) — 0.80. The structural con- 
figurations associated with the crystallization transition, 
as a function of filling fraction, were studied in detail 
in and the caging dynamics, as crystallization is ap- 




FIG. 3: Snapshots of typical configurations of the granular 
layer at various values of the filling fraction: (a) dilute gas, 
(j> = 0.34, (b) dense liquid, = 0.67 and (c) thermalized 
crystal, = 0.80. Driving parameters: / = 50Hz, T = 4.0. 
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FIG. 4: Filling fraction dependence of the granular temper- 
ature {Tx and Ty are the x and y components, respectively), 
for both cases of using a flat and rough bottom glass plates. 
Driving parameters: / = 50Hz and F = 4.0. 



proachcd in Q. 

As defined in the previous section the granular temper- 
ature is the average kinetic energy per particle. For this 
the brackets (.) in Eqn. (|16p now denote both time and 
ensemble averages for all the spheres found within the 
imaging window. Moreover, one can define — l/2(u^) 
and Ty = l/2{vy), with T — T^+Ty, as the temperature 
projections onto the x and y directions of the horizontal 
plane. 

In Fig. |4]the temperature components, T^ and Ty, are 
plotted as a function of filling fraction. For both cases 
of using a rough and a flat bottom plate, the respective 
Tx and Ty are identical. This shows that the forcing of 
the granular particles was isotropic. In agreement with 
the case of a single particle, the dynamics of the granular 
layer using a flat or a rough bottom plate is remarkably 
different. If the rough bottom plate is used, the gran- 
ular temperature depends almost monotonically on the 
filling fraction. At low cj), T is approximately constant 
as the layer simply feels the thermal bath and shows lit- 
tle increase in T until (p ~ 0.5 is reached. As the filling 
fraction is increased past (j) ~ 0.5, the granular tempera- 
ture rapidly decreases due to energy loss in the increas- 
ing number of collisions and due to a decreasing available 
volume. 

For the case of using a fiat bottom plate, the non- 
monotonic dependence of T on is more dramatic and it 
is difficult to attain homogeneous states below (j) < 0.4. 
This is due to the fact that at those filling fractions 
(the limiting case being the single particle investigated in 
the previous section) the particles perform small excur- 
sions and interact with their neighbors only sporadically. 
Hence, little momentum is transferred onto the horizon- 
tal plane. For (j) > 0.5, there is an increase of temper- 
ature with increasing filling fraction as interaction be- 
tween neighbors becomes increasingly more common up 



to ~ 0.7 at which the curve for the flat bottom plate 
coincides with the curve for the rough bottom plate. It 
is interesting to note that the value at which this match- 
ing occurs is close to t he p oint of crystallization of disks 
in 2D, = 0.716 [1 SJ Sll. At this point, the large 
number of collisions between neighboring particles ther- 
malizes the particles, independently of the details of the 
heating, i.e. of whether a flat or rough plate is used. 

The behavior of the layer was therefore more robust 
by using a rough bottom plate. Moreover, this allowed 
for a wider range of filling fractions to be explored, in 
particular in the low limit. Further detailed evidence 
for the advantage of using the rough bottom plate over 
a flat one is given in Appendix [A) Hence, all results 
presented in the remainder of this paper correspond to 
experiments for which a rough glass plate was used for 
the bottom plate of the experimental cell. 

In Fig. [5] we present the dependence of the granular 
temperature of the layer on the forcing parameters for 
a fixed value of filling fraction, (j) ~ 0.59. There is a 
monotonic increase of T with both / and F. 

The final control parameter that we investigated was 
the height of the experimental cell which we varied by 
changing the thickness of the inter-plate annulus by us- 
ing precision spacers. The range explored was 1.3D < 
h < 2.3D. Note that for h ^ D there would be no clear- 
ance between the spheres and the glass plates and one 
would therefore expect for no energy to be injected into 
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FIG. 5: Dependence of the granular temperature, T, on driv- 
ing parameters: (a) frequency, / and (b) dimensionless accel- 
eration F. The filling fraction is kept constant at (/> = 0.59. 
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FIG. 6: Dependence of the granular temperature, T, on cell 
gap, h, which was varied by changing the thickness of the 
inter-plate annulus. Filling fraction was kept constant at 
o!)=0.59 and the dimensionless acceleration F was fixed to 4. 



the system. On the other hand for h > 2D, spheres 
could overlap over leach other and the system ceases to 
be quasi-2D. The dependence of the granular tempera- 
ture on cell gap is plotted in Fig. [6l It is interesting to 
note that T appears to depend linearly on the gap height. 

One can therefore regard changing h, f and F as a way 
of varying the temperature of the granular fluid, which 
can this way be tuned up to a factor of eight. 



VI. PROBABILITY DENSITY FUNCTIONS OF 
VELOCITIES 

We now turn to the distribution functions of particle 
velocities under various conditions of filling fraction, fre- 
quency, dimensionless acceleration and gap height. In 




FIG. 7: (a) Probability distribution function of velocities 
(PDF), P{v), for specific values of filling fraction {(p), fre- 
quency (/) and dimensionless acceleration (F). See legend for 
specific values. 



Fig. [7] we plot the probability density function of a ve- 
locity component, P{vi), where the index i represents 
the component x or y, for specific values of 0, / and F. 
From now on we drop the index in Vi since we showed in 
the previous section that the dynamics of the system is 
isotropic in x and y. Moreover, when we refer to velocity 
we shall mean velocity component, i. The width of the 
observed P{v) differs for various values of the parame- 
ters. This is analogous to the fact that the temperature 
(i.e. the variance of the distribution) is different for vari- 
ous states set by <f), f and F, as discussed in the previous 
Sec. It is remarkable, however, that all the P{v) distri- 
butions can be collapsed if the velocities are normalized 
by the characteristic velocity, 



(17) 



where the brackets (.) represent time and ensemble av- 
eraging over all the particles in the field of view of the 
imaging window. This collapse is shown in Figs. [H^a) 
and (b) and was accomplished not only for various filling 
fractions but also for a range of frequencies and dimen- 
sionless accelerations. 

To highlight the quality of the collapse we have sep- 
arately plotted P{c) for the four lowest values of filling 
fraction (0 = 0.27,0.32,0.34,0.39) along with the P(c) 
for a single particle in Fig. [S^a) and all other data in 
the ranges 0.44 < < 0.8, AOHz < f < lOOHz and 
1.4 < F < 6 in Fig. [S^b). At low values of the filling 
fraction (0 < 0.44) the collapse is satisfactory but devi- 
ations arc seen at low c. In particular, near c — the 
distributions exhibit a sharp peak with a clearly discon- 
tinuous first derivative which reflects the fact that for 
these value of low (j) the resulting gases are not collision- 
ally driven but are dominated instead by the underlying 
thermostat (c.f. distribution for single particle). How- 
ever, this sharp peak becomes increasingly smoother as 
the filling fraction is increased, presumably due to the 
increasing number of particle collisions, and by ^ 0.44, 
it has practically disappeared. On the other hand, for 
4> > 0.44 the collapse onto a universal curve is remark- 
able for such a wide range of 0, / and F. The choice of the 
lower bounds for frequency and acceleration, / = AOHz 
and F = 1.4, in the data plotted will be addressed in Sec. 
IVIBI Moreover, this collapse of P{v) is also attained for 
various values of the gap height {1.3D < h < 2.3D), as 
shown in Fig. ^c). From now onwards, we shall perform 
our analysis in terms of the reduced velocity, c = v/vq- 

We stress the universal collapse of the experimental ve- 
locity distribution functions for a wide range of parame- 
ters and the clear deviation from the Maxwell-Boltzmann 
(Gaussian) distributions of Eqn. (fTO| - solid curve in Fig. 
^h) and (c). This is particularly visible at high velocities 
where there is a significant overpopulation of the distri- 



bution tails in agreement with previous theoretical 33 1 , 



numerical [3^ and other experimental work [27|, [3] 
To quantify these deviations from Gaussian behavior, our 
analysis will be twofold. First we shall analyze the non- 
Gaussian tails of the velocity distributions (Sec. IVI Ap . 
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Even though the tails correspond to events associated 
with large velocity the probability of them happening is 
extremely low. However, this overpopulation at the tails 
implies that, if the distribution is to remain normalized 
with a standard deviation given by the average kinetic en- 
ergy per particle, clear deviations from Gaussian at the 
central regions of low c (high probability) of the distribu- 
tion must necessarily also be observed. In Sec. IVIBI we 
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FIG. 8: (a)-(b)PDF of velocities, -P(c), in which all velocities 
for each distribution were normalized by its standard devia- 
tion. Data for a range of (p, f and F. (c) PDF of velocities 
for range of values of cell gap, h. Other parameters set to 
(j> = 0.59, / — 50Hz and F = 4.0. (a-c) See legend for specific 
values of the experimental parameters. The solid lines are 
Gaussian with unit standard deviation. 



study this deviations in the context of Sonine corrections 
introduced in Sec. |TT1 

In Appendix we provide the data corresponding to 
Fig. [7] and Figs. [Ha,b) of the velocity distributions ob- 
tained using the optically flat bottom plate, instead of the 
rough plate used here. We shall show that the P{v)s for 
the flat bottom plate exhibit much stronger non-universal 
deviations from Gaussian behavior and that the collapse 
of P{c) is highly unsatisfactory. This highlights, as men- 
tioned in Sec. |Vl the considerable advantage in our ex- 
perimental technique of using the rough plate to generate 
the granular fluid. 



A. Deviation at large velocities — The tails 

Recently, van Noije and Ernst 3^ have made the the- 
oretical prediction that the high energy tails of the veloc- 
ity distributions of a granular gas heated by a stochastic 
thermostat should scale as stretched exponentials of the 
form, 



Pic) 



exp 



(18) 



where A is a fitting parameter. Note that the theoret- 
ical argument of van Noije and Ernst involves a high c 
limit and hence one does not expect this stretched ex- 
ponential form to be valid across the whole range of the 
distribution. 

To check the applicability of Eqn. (fTS]) to describe 
our data, we have calculated the quantity q(c) = 
— In (— ln(P(c))), which we plot in Fig. ^a,) for a va- 
riety of filling fractions, frequencies and dimensionless 
accelerations. Indeed, within the ranges considered and, 
in the limit of large c, g(c) tends to a straight line with 
slope -3/2 for all values of the control parameters, in ex- 
cellent agreement with the scaling of Eqn. (fT8| . In Fig. 
mja) we have excluded filling fractions < 0.44, since 
as discussed in the previous Sec, the gases obtained in 
this range are dominated by the thermostat rather than 
collisions and the main ingredient for the prediction of 
Eqn. (fT5|) is the role of inelastic collisions. This behavior 
of the tails of the velocity distribution is also observed 
for various values of the gap height, h, of the experimen- 
tal cell, as shows in Fig. Hfb). Note the crossover from 
Gaussian-like (dashed line with slope -2) to the stretched 
exponential (solid line with slope -3/2) at c ~ 1, which is 
particularly clear at the larger values of 0. For low (j> the 
tails of the distribution tend to be closer to the stretched 
exponential throughout the range —1.1 < log(c) < 1.1. 
We stress that Eqn. (flS)) has been derived in the limit of 
large c and is therefore only expected to be valid at the 
large c end of the tails. A. Barrat el. al. [44, 45] have 
investigated the range of validity in a numerical model 
with a stochastic forcing and found that the stretched 
exponential scaling of the tails should only be valid for 
very high velocities corresponding to probabilities lower 
than 10~^P(0). The fact that we find this scaling at con- 
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/a/s, such that, 



FIG. 9: (a) Tails of P{c) for specific values of filling frac- 
tion {(f}), frequency (/) and dimensionless acceleration (F) and 
(b) gap height (ft). The solid lines correspond to stretched ex- 
ponentials of the form ~ exp{—Ac'^^^) whereas the dashed 
lines correspond to the Gaussian behavior of the form ~ 
exp(-Ac^). 



siderable lower velocities, i.e. in regions of higher proba- 
bility, (in agreement with other experimental work [2^) 
raises some questions regarding the detailed comparison 
between theory and experiments. 



B. Sonine corrections to the distribution 

Having looked at the tails of the distributions, i.e. at 
large c, we now extend our deviation analysis based on an 
expansion method to the full c range of the distributions. 
This highlights in particular deviations from Gaussian 
behavior in the central high probability regions of P{c) 
near c = 0. In Sec. [IT] we discussed that in the solution 
of the Enskog-Boltzmann equation for inelastic particles 
in a stochastic thermostat, a Sonine expansion is usu- 
ally performed such that the deviations from Gaussian 
are described by a Sonine Polynomial (i.e. a 4th order 
polynomial with well defined coefficients) multiplied by 
a numerical coefficient a2- 

To check the validity of this assumption of the Kinetic 
Theory, we shall follow a procedure analogous to that em- 
ployed in the numerical study of Ref. [S^j . We calculate 
the deviation, A(c), of the experimental velocity distri- 
butions, -P(c), from the equilibrium Maxwell-Boltzmann, 



P{c) = fMB (1 + A(c)) 



(19) 



where Jmb is given by Eqn. pO|) . By comparing Eqn. 
(fin)) to Eqn. ^ for the theoretical case of inelastic par- 
ticles under a stochastic thermostat, we expect that 



i"(c) = /a/s(1 + 02^2(0')), 



(20) 



such that the experimental deviations from equilibrium 
take the form of the Sonine polynomial of order-two, 
S2{c'), i.e. 



A(c) = 02 (l/2c^ - 3/2c2 + 3/8) 



(21) 



where the parameter 02 can be directly related to the 
kurtosis K of the distribution as 02 = K/3 — 1. However, 
because of sampling noise at high c, we have taken 02 as a 
fitting parameter rather than determining it directly from 
K, although in most cases there is little difference. Note 
that the fitting of the experimental data to Eqn. pp]) 
was done over the whole range of the distribution with a 
relative weighting given by the corresponding probabil- 
ity. In Fig. [TUTa-d) we plot these experimental deviation 
from equilibrium, A(c), for four representative values of 
filling fraction, along with the order-2 Sonine polynomial 
(solid curve) as given by Eqn. (PT|) and fitting for a2. At 
low filling fractions, for example cj) = 0.34 - Fig. [TCTa) - 
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FIG. 10: Experimental deviation function from Gaussian be- 
havior, A(c) for (a) = 0.34, (b) = 0.53, (c) cj> = 0.64 and 
(d) (j) — 0.80. The solid line is the order-two Sonine polyno- 
mial of the form a2(l/2c'' — 3/2c^ +3/8) where 02, the second 
Sonine coefficient, is the only adjustable parmeter . 
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FIG. 11: Experimental deviation function from Gaussian be- 
havior, A(c) for (j> — 0.66. The sohd hne is the Sonine polyno- 
mial 02 (l/2c'* — 3/2c^ -1-3/8) with one single fitting parameter: 
a2 = 0.171. The dashed line is the higher order Sonine poly- 
nomial description of the form Ylp=2 <^pSp{(?) with the follow- 
ing (five fitting parameters) Sonine coefficients; 02 = 0.1578, 
as = -0.0656, 04 = 0.1934, as = -0.1637 and ae = 0.0832. 



the second order Sonine fit is approximate but unsatis- 
factory for small c where a sharp cusped peak is present. 
This cusp is reminiscent of the underlying dynamics of 
the rough plate thermal bath as evidenced by the cusp 
at c = for the single particle distribution presented in 
Fig- EIc) and discussed previously. On the other hand, 
for larger values of the filling fraction (in particular for 
<j) > 0.44) the experimental data is accurately described 
by the second order Sonine polynomial. 

To further evaluate the relevance of the Sonine poly- 
nomials to describe the experimental deviation from 
Maxwell-Boltzmann, in Fig. [TT] we plot the experimen- 
tal A(c) along with both the order-two {p = 2) Sonine 
polynomial and the higher order Sonine expansion of the 
form X]p=2 '^p'^p('^^)- second order Sonine polyno- 

mial alone is responsible for the largest gain in accuracy 
of the corrections from MB (horizontal dashed line). The 
higher order Sonine expansion indeed provide a better fit 
but the overall improvement is only marginal (see Fig. 

El) 

In Fig. [T2ja) we present the dependence of the sec- 
ond Sonine coefficient a2, the single fitting parameter in 
Eqn. (|2T|) . as a function of filling fraction. The coef- 
ficient a2 initially decreases with increasing filling frac- 
tion wp to (f> — 0.65, after which it shows a rapid rise. 
We also study the dependence of a2 on the driving fre- 
quency (/) and dimensionless acceleration (T) for a single 
value of the filling fraction (p = 0.59, which we present 
in Figs. [iWb) and (c), respectively. For lower values 
of frequency (20-40Hz) we get a monotonic drop in 02 
with increasing frequency. For / > AOHz, the coefficient 
a2 then levels off and remains approximately constant 
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FIG. 12: Experimentally determined order-two Sonine coef- 
ficient a2 as a function of (a) filling fraction (•), (b) driv- 
ing frequency (□), (c) dimensionless acceleration (A) and (d) 
granular temperature. The filling fraction is kept constant at 
(j) = 0.59 while exploring the dependence of 02 on / and P. 



at 02 0.171 ± 0.023. For F > 2.0, the coefficient 02 
remains approximately constant at 02 — 0.193 ± 0.029. 
We have plotted these two datasets at fixed (p and vary- 
ing / and r back in Fig. [TWa) and the behavior of 02 
is found to be consistent with the previous dataset with 
varying cj) at fixed (/, F) = (SO-ffx, 4.0). It is interest- 
ing to note that a scatter of points is obtained if 02 for 
the three previous datasets is plotted as a function of 
the corresponding granular temperatures. Moreover, for 
the data-set for varying (j) the dependence is clearly not 
single valued, i.e. one can find two distinct values of 02 
for a single value of temperature. All this suggests that 
the temperature does not set 02. Instead, 02 is a strong 
function of (j) and is only weakly dependent on / and F, 
provided that / > AOHz and F > 2. This is in agreement 
with inelastic hard-sphere behavior where the only two 
relevant parameters are thought to be the filling fraction 
and the coefficient of restitution (not investigated in the 
present study). 

Next, in Fig. [T3] we present the dependence on fill- 
ing fraction of the experimentally determined five non- 
zero Sonine coefficients in the sixth order expansion of 
^(c) ^ J2p=2'^p'^pi'^'^)- III theoretical analysis, the 
higher order terms (order 3 and above) in the Sonine 
polynomial expansion are typically neglected to simplify 
the calculations, the claim being that this results in no 
significant loss of accuracy. In our experimentally gener- 
ated fluids the higher order Sonine coefficients assume 
values less than 1 for flUing fractions above = 0.5. 
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FIG. 13: Dependence of higher order coefficients of the Sonine 
expansion (02 to ae) on filhng fraction. 



From this it is clear that it is unnecessary to consider 
orders higher than two for intermediate and high fiUing 
fractions. For (j) < 0.5, however, they are significant and 
the finite values of 03,..., ae are required to fit the in- 
creasingly larger regions of stretching exponential tails 
which become progressively more wider and propagate 
towards smaller velocities. As we mentioned in the dis- 
cussion of Fig. [Ifa) in the previous Sec, when the filling 
fraction decreases, the dynamics of the layer starts re- 
sembling that of the underlying thermal bath set by the 
rough plate. There the distribution, at low c becomes 
increasingly more like a stretched exponential through- 
out the full c range, rather than only in its tails as is the 
case at larger values of (j). This finding is analogous to a 
result from Brilliantov and Poschel ^32] who found both 
theoretically and numerically that the magnitude of the 
higher-order Sonine coefficients can grow due to the an 
increasing impact of the overpopulated high-energy tails 
of the distribution function, just like in our system at 
low (f). Their analysis was however performed as a func- 
tion of the coefficient of restitution, e, rather than fill- 
ing fraction. In their study a single second order Sonine 
polynomial becomes insufficient to represent the devia- 
tion from Gaussian at low values of e. Moreover, they 
worked on the Homogeneous Cooling State (HCS), but 
they suggest that their results should also apply to gran- 
ular gases with a thermostat, like in our case. A more 
direct comparison with this theoretical analysis is open 
to further investigation. 

We stress that a finite value of the higher order coef- 
ficients docs not necessarily imply a large correction to 
the deviations A(c). To explore further this point we 
now quantify the deviations of our experimental velocity 
distribution from all the three models we have consid- 
ered: 1) the equilibrium Maxwell-Boltzmann distribu- 
tion, 2) the velocity distribution function with order-2 
Sonine polynomial expansion and, 3) the velocity dis- 
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FIG. 14: Deviations (quantified by x^)i a function of fill- 
ing fraction, of the experimental velocity distributions from 
Maxwell-Boltzmann distribution (□), from the velocity dis- 
tribution function with a order-two Sonine polynomial cor- 
rection (v) and from the velocity distribution function with 
a order-six Sonine polynomial expansion correction (o) . 



tribution function with higher order Sonine polynomial 
terms. In Fig. [T3]we plot these deviation of experiments 
from the models as quantified by = J2iLiiPi^i)exp — 
P{ci)modei)'^ , for the full range of filhng fractions. We 
clearly see significant deviation of the experimental data 
from Maxwell-Boltzmann distribution for all (j), whereas 
the velocity distributions with the order-2 Sonine poly- 
nomial term shows a considerably better agreement with 
the experimental data across the whole range of (j>. Even 
though the higher order Sonine polynomial expansion 
characterizes the experimental data more closely, these 
higher order contributions are modest compared with the 
large gain in accuracy from the order-2 Sonine polyno- 
mial alone, and in the range of (/) > 0.44 provide no sig- 
nificant improvement. Thus the order-2 description ob- 
tained by neglecting the higher order Sonine terms is a 
reasonably good approximation to describe the experi- 
mental data. 



VII. CONCLUSION 

In conclusion, we have developed an experimental 
model system for a quasi-2D granular layer, under homo- 
geneous stochastic driving. Our experimental technique 
has allowed us to randomly thermalize a granular fiuid 
over a wide range of filling fractions. Our study was cen- 
tered on the dynamics of the system, in particular the 
statistics of velocities. The temperature of the experi- 
mental granular fluid could be adjusted by varying the 
system's control parameters, namely the filling fraction, 
the frequency and acceleration of the driving and the gap 
height, with no significant change in the nature of the ve- 
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locity probability distribution functions. We have found 
an excellent collapse of the distribution functions if the 
particle velocities are scaled by a characteristic velocity 
Vo (the standard deviation of the distributions). 

However, the obtained distribution are non-Gaussian. 
We have analyzed the deviations from Gaussian behavior 
in two distinct regimes. Firstly, we looked at the shape of 
the distribution tails which scaled as stretched exponen- 
tials with exponent —3/2. Secondly, we have performed 
an expansion method that highlights the deviations from 
a Gaussian at low velocities (near c — 0) and found them 
to be Sonine-like, i.e. polynomial of order four with fixed 
coefficients. This way we have determined the validity 
of some important assumptions in the Kinetic Theory 
of randomly driven fluids. It is surprising that this for- 
malism seems applicable at filling fractions as high as 
(j> — 0.80, whereas Kinetic Theory is usually thought to 
breakdown at much lower values of (f). Finally, we have 
looked at Sonine polynomial expansion with higher order 
terms and concluded that it is sufficient to retain only the 
leading order (first non-vanishing) term to maintain rea- 
sonable accuracy of the polynomial expansion. There- 
fore, we can accurately characterize the single particle 
velocity distribution function by introducing a single ex- 
tra coefficient 02, in addition to the 0th, 1st and 2nd 
moments used for fluids at equilibrium. The coefficient 
02 has universal character in the sense that it is a strong 
function of filling fraction and only depends weakly on 
the other experimental parameters such as the driving 
frequency and acceleration. 

Having determined the base state of our randomly 
driven granular fluid, we hope that an experimental sys- 
tem such as ours can be used further as a laboratory 
test bed of some basic assumptions of Kinetic Theory. 
In future work, it would also be of interest to investi- 
gate whether particle-particle velocity correlations [2§| 
are present in our system but this is beyond the scope of 
the current investigation. 

To our knowledge, this is the first time that the So- 
nine corrections of the central high probability regions 
of the velocity distributions have been measured in an 
experimental system, in agreement with analytical pre- 
dictions. This should open way to further theoretical 
developments, crucial if we are to develop much desired 
predictive models for granular flows with practical rele- 
vance. 



APPENDIX A: PROBABILITY DENSITY 
FUNCTIONS OF VELOCITIES USING THE 
BOTTOM FLAT PLATE 

In this Appendix we discuss the nature of the veloc- 
ity distributions obtained using a flat bottom plate as in 
the system of Olafsen and Urbach [2^ . The distributions 
P{v) are presented in Figs. [TSTa) for a wide range of flU- 
ing fractions. Unlike the case of using a rough plate, here 
there is an enormous variation in the shape of Piv). At 



low 4> the distributions are sharply peaked near c = 0. 
The value of P{v — 0) can differ by over a factor of 20 
between low and high (f>. In Fig. [THb) the distribu- 
tions P{c) - obtained after rescaling all velocities by the 
characteristic velocity Vo given by Eqn. p7|) - do not 
collapse onto a single curve and show considerably larger 
deviations from a Maxwell-Boltzmann (solid line). This 
distributions only attain the general form found for the 
rough case at filling fractions larger than (j> > 0.6 and this 
is in agreement with the point at which the granular tem- 
perature of both cases coincide - Fig. [H For this high 
values of (j>, both in the rough and flat cases, the con- 
tribution to the thermalization is more important from 
the large number inter-particle collisions rather than the 
underlying thermostat. Hence we expect the analysis we 
have performed throughout this paper to be applicable 
to the flat case, for (f> > 0.6, but not below this point. 

In addition to the discussion on the granular tempera- 
ture of Sec. |Vl these results present further evidence for 
the significant advantage of using a rough bottom plate 
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FIG. 15: PDF of velocities for an experimental cell where 
both top and bottom plates were optically flat, (a) P{v) (b) 
P{c), in which all velocities for each distribution were normal- 
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Boltzmann distribution with unit standard deviation. Data 
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to generate the thermostat to heat the granular fluid in 
a way that may be directly compared to theoretical work 
on inelastic hard-spheres driven by a stochastic thermo- 
stat. 
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